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In the present work, we consider the self-focusing discrete nonlinear Schrodinger equation on 
hexagonal and honeycomb lattice geometries. Our emphasis is on the study of the effects of 
anisotropy, motivated by the tunability afforded in recent optical and atomic physics experiments. 

We find that important classes of solutions, such as the so-called discrete vortices, undergo desta¬ 
bilizing bifurcations as the relevant anisotropy control parameter is varied. We quantify these 
bifurcations by means of explicit analytical calculations of the solutions, as well as of their spectral 
linearization eigenvalues. Finally, we corroborate the relevant stability picture through direct nu¬ 
merical computations. In the latter, we observe the prototypical manifestation of these instabilities 
to be the spontaneous rearrangement of the solution, for larger values of the coupling, into local¬ 
ized waveforms typically centered over fewer sites than the original unstable structure. For weak 
coupling, the instability appears to result in a robust breathing of the relevant waveforms. 


I. INTRODUCTION 

In both optical media [lj] and atomic systems, such as Bose-Einstein condensates (BECs) [2i] , in the past two decades 
there has been a tremendous amount of effort focused on understanding the implications of periodic lattices. In the 
former case, both the realms of optical waveguides [|] and of photorefractive crystals Q have played crucial roles 
towards the analysis and experimental realization of states such as discrete solitons, and vortices, as well as of more 
complex waveforms, including ring structures, necklaces, gap solitons and many others. In atomic BECs, on the other 
hand, the emphasis has not only been on corresponding matter waves jf|, but also on quantum phenomena beyond 
the realm of mean-field models [6j. 

In recent years, the emphasis has somewhat shifted from the consideration of the more customary square lattices to 
the examination of lattices of hexagonal or honeycomb form. There, a source of emphasis has again been localization 
and self-trapping in the form of solitonic and vortical structures but also other aspects have been studied 

including, e.g., Bloch states 0 - A significant fraction of the focus has been on the emulation by these optical systems 
of “photonic graphene”, leading to numerous remarkable features^including the creation, destruction and experimental 
observation of topologically protected, so-called, edge states [lllllMl, and also the emergence of pseudospin and angular 
momentum [13]. In the atomic realm too, considerably tunable and flexible optical lattices of both a hexagonal and 
honeycomb form have been produced for single [l4| and multi-species [15] experiments. While much of the interest 
in this context lies within quantum mechanical transitions, such as the superfluid-insulator transition 16 ], the atoms 
can, very controllably, be considered in the superfluid regime where a mean-field description paralleling the optical 
one is suitable. As an aside, it is relevant to mention that more complex lattice structures including e.g. Kagome 
lattices are also a subject of ongoing consideration 00 and are within the realm of experimental possibility in 
both settings. 

At the mathematical level, there exists a prototypical model that combines the suitable lattice geometry, the 
discreteness and the nonlinearity. As a result, it captures the principal features of the experimental observations, at 
least as regards the emerging coherent structures. This model is the so-called discrete nonlinear Schrodinger (DNLS) 
equation, which has been a subject of intense theoretical and numerical investigation [ 1 8l ] . Our aim in the present 
work is to utilize this DNLS model in order to capture the impact of anisotropy on the hexagonal and honeycomb 
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lattices. This is in part motivated by the studies in optical photorefractive systems such as the work of [§] where 
both unstretched and stretched lattices were used and in both cases the coupling was anisotropic (varying in one 
direction between 20% and 80% of the coupling in the other directions). Such a systematic study is also motivated 
by the atomic realm of, e.g., [Tij . where the full control of the optical beam intensities, wavenumbers and phases that 
create the lattice trapping the atoms can straightforwardly be used to produce different types of lattices (e.g. both 
hexagonal and honeycomb) and different anisotropies. 

Our aim here is to provide a systematic analysis of the different types of solutions that are possible in the anisotropic 
system. Starting from the isotropic two-dimensional limit, we vary the strength of the interaction along a particular 
direction. Progressively this leads from a two-dimensional configuration, e.g. in the honeycomb case, to an uncoupled 
set of quasi-one-dimensional configurations. As a result, we can appreciate that numerous states among those that 
exist in the two-dimensional (2d) setting should disappear at a suitable critical point as we approach the Id regime. 
For instance, the discrete vortices belong to this category, as there are no solutions with nontrivial vorticity in one¬ 
dimensional DNLS lattices mm. Here, we intend to provide a quantification of the relevant solutions, as well as to 
provide a road map for their dynamical destabilization by evaluating their dominant linearization eigenvalues. Both 
of these steps are performed analytically (to leading order) permitting a complete characterization of the bifurcation 
events/destabilization or disappearance of different branches of solutions. This is done for the prototypical unit cell of 
each lattice i.e., for a triangular cell within the hexagonal lattice and a hexagonal cell within the honeycomb lattice, 
although it can be straightforwardly generalized to other cases. Once the existence, stability and bifurcations of 
the relevant solutions are determined, then their potential instabilities (and spectral properties) are also explored 
numerically. Finally, these findings are corroborated by direct numerical computations illustrating the tendency of 
the (unstable) dynamics towards (typically) fewer sites than the original structure. In the case of the 3-site cell in 
the hexagonal lattice, we observe a tendency of the dynamics towards the ground (single-site) state of the model for 
stronger couplings, or towards robust breathing excitations in the case of weaker couplings. In the case of the 6-site 
cell of the honeycomb lattice, even for stronger couplings, multi-site excitations (of different types - see details below) 
were typically found to persist over the evolution scales of dynamical propagation considered here. 

Our presentation is structured as follows. In section II, we present our systematic analytical findings regarding 
the existence and stability of solutions for each of the lattices in their respective unit cells. Then in section III, we 
present the corresponding numerical findings, as well as examine, the fate of dynamically unstable solutions. Finally, 
in section IV, we summarize our results and present some challenges for future studies. 


II. THEORETICAL ANALYSIS 


To study the two geometries of interest, we consider the following discrete nonlinear Schrodinger equation 


.du„ 


— = -£A- 2 Um,n - \u. 


dz 

with the two-dimensional anisotropic discrete Laplacian 


m.n I u j m,n 
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( 2 ) 


Here, the constant e denotes the strength of the linear coupling between nearest neighbor sites in the isotropic case 
and the anisotropy is controlled by the parameter 6 (0 < S < 1). Since both of the grids of interest are rotationally 
invariant, hence there is a freedom in selecting the direction of the anisotropy. A value of S = 1 yields the isotropic 
lattice with a uniform nearest neighbor coupling of e, whereas 5 = 0 completely decouples sites in the direction parallel 
to a particular lattice direction, chosen without loss of generality. Physically, the field represents (in our optical 
example) the envelope of the electric field in the waveguide, while in that case z is the propagation coordinate. For 
BECs, the field represents the atomic wavefunction in the corresponding well of the optical lattice, while z in that 
realization is replaced by the time t. The summations are over disjoint subsets, Ni and N 2 , of the set N = NiU N 2 
of nearest neighbors where |AT| = 6 for the hexagonal lattice and \N\ = 3 for the honeycomb lattice. The set N 2 is 
the set of nearest neighbors joined to u m ^ n by the (anisotropic) coupling, Se , while the remaining nearest neighbors 
belong to the set N 2 and have a coupling of e to Note that in the case 5 = 0, the hexagonal grid becomes 

the usual rectangular grid, while the honeycomb grid becomes a parallel set of one-dimensional grids with (in both 
geometries) a nearest neighbor coupling of e. 

We are interested in stationary solutions of the form u m ^ n = exp(iAz)ri m „, where A is the propagation constant in 
optics or the chemical potential in BECs. Then, v m , n satisfies the steady-state equation 

AVm^n = ~\~ |^r?T,, 7 T,| ^771,n- 


( 3 ) 
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In the anticontinuum (AC) limit of uncoupled sites [2(j, (i- e - when e —> 0), the solutions of Eq. ([3]) are u mjra = 0 
and v m , n = -\/Aexp(i0 mjn ). Thus, at the AC limit, explicit solutions of the form u k = -\/Aexp(i0*.) exp(iAz) can be 
found over contours M of the uncoupled lattice points for arbitrary 9 k E [0, 27 t), where the nodes are indexed by k. 
For the current work, in the hexagonal lattice, k will index the sites along a three-site one-dimensional closed contour 
(\M\ = 3), while in the honeycomb lattice, k will index the sites along a six-site one-dimensional closed contour 
(\M\ = 6), i.e., we will consider the principal cells of the respective lattices. Without loss of generality, we set A = 1. 

Following an analysis similar to that of m, mi, we then find that the necessary (leading order) conditions for 
solutions over a discrete contour to persist for e > 0 are given by 

F k = 5k,k -1 sin (8 k - 9 k -i) + 5 k , k +i sin(0 fc - 9 k+1 ) = 0, (4) 

where we have the periodic condition 9 k+ \M\ = 6k for k = 1,..., \M\, and the coefficients 5 k , k -i and 5 k , k +i, provide 
the lattice anisotropy as defined by: 


5 k ,i 


5, if the segment of M joining adjacent nodes k and l is parallel to the anisotropic direction 
1, if the segment of M joining adjacent nodes k and l is not parallel to the anisotropic direction. 


(5) 


We will study the behavior of some solutions of the variety described above in both the hexagonal and honeycomb 
lattices when the “background” coupling, e, is fixed at a small value and the anisotropy is switched on (i.e. 6 < 1 ). 
We will see that the anisotropy may stabilize or destabilize some solutions via, typically, pitchfork (i.e., symmetry 
breaking) bifurcations. Theoretically we will find that in the weak coupling limit, this transition from stability to 
instability, or vice-versa, occurs when the solution collides with another solution at <5 = 0.5 which then persists as the 
lattice becomes more anisotropic, until there is a complete decoupling in the prescribed direction. 

Again, adapting the results of mm (see also the exposition of }18j). the stability of lattice excitations can be 
determined (to leading order) from the eigenvalues j k of the \M\ x \M\ Jacobian matrix of the form 


{ 5 k , k —i cos {9 k 0fc_i) “t - cos (9 k 9 k -\. i), k — l 

—6 k ,icos(6 k — 6i), k = l± 1 ( 6 ) 

0 , \l — k\> 2 . 

When the excited nodes in the lattice are adjacent, the (near zero) stability eigenvalues X k of the full problem are 
given by X k = ±^/2') k e [H, l2lj |. We now examine some explicit examples of this general theoretical formulation of 
the anisotropic DNLS problem. 


A. Hexagonal Lattice 


In this subsection, we will analytically track the effects of anisotropy on the stability of various three-site configu¬ 
rations in the hexagonal lattice and make predictions about what the original configurations in the isotropic lattice 
transform into. The notation [a, b , c] is employed to describe the three-site contour with phases a, b and c at the 
nodes (while a corresponding notation will be used for six-site contours). The contours we discuss here are: (1) 
[0, 27 t/3, 47 t/ 3 ] (charge one vortex), (2) [0, n, 0], and (3) [0, 0,0]. These are the principal (up to trivial transformations 
of phase) 3-site solutions of the hexagonal lattice cell. 

(1) We begin with the three-site single charged vortex (0i = 0, |A0| = ^). From Eq. Q with \M\ = 3 and with 
the anisotropy lying between the sites with phases 8\ and 0 3 , we obtain the following relationship between the phases: 


9 3 = 9 1 + 2 arccos ( - — 


+ i 


do = 


Using this phase profile in Eq. ([!>]), the Jacobian matrix becomes 


J = 


'0i-h. 


— COs(- 2 J 

—5 cos (0i — 9 3 ) 


0.5 < <5 < 1 
(mod 2n). 


cos ( 9l 2 e3 ) + ^cos(0i — 0 3 ) — cos ( 01 2 03 ) 

2cos(^) 


—6 cos ( 0 i — 03 ) 
-cos(^) 


- COS ( 01 2 03 ) COS ( 01 2 6,3 ) + J COS (01 - 0 3 ) 


(7) 


( 8 ) 
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Thus, for 0.5 < 5 < 1 the eigenvalues 7 fc are found to be 71 = 0, 72 = cos ( 81 2 03 ) + 25cos($i — 63) and, 73 = 
3 cos ( 9l ^ 93 ). In particular, for the hexagonal three-site charge-one vortex (which satifies Eq. (|7|)) the stability 

eigenvalues for 0.5 < 5 < 1 are found to be Ai = {0,0}, A 2 = +\j and A 3 = For this range of 

values of 5, all of the quantities in the radicals within the eigenvalue pairs remain nonnegative and, for sufficiently 
small e, the pair A 3 will not collide with the continuous spectrum. Hence, this vortex remains stable throughout this 
interval of anisotropy. Note that the only eigenvalue pair that moves along the imaginary axis towards the origin of 

the spectral plane and thus has the potential to bring about instability when 0 < <5 < 0.5 is A 2 = ± \J As 

the above solution of Eq. m cannot be continued below 5 = 0.5, additional analysis is needed to reveal the outcome 
for a further increase in anisotropy, i.e. for 0 < 5 < 0.5. 

From Eq. 0, we can also determine the theoretical predictions for the changes in the relative phases (mod 27 t) as 
a function of 5 for 0.5 < 5 < 1: 


\d 2 - 6 x\ = \9 3 - 0 2 \ = ^| 6 »i - 0 3 1 = arccos • ( 9 ) 

In the isotropic case (i.e. when 5 = 1) for all the relative phases we have |A(9| = As the anisotropy increases, 
reaching 5 —► 0.5, we find that \0 2 — 6 *i| —>■ 7 r, \d 3 — d 2 \ —» n, and \9\ — 0 3 \ —» 0. Thus, at this critical point, the 
single-charge vortex [ 0 , ^L], merges with the configuration [ 61,61 + 73,63 + 2 n\ = [ 61,63 + n, 0 i](mod 27r) i.e., with 

the [0, 7 r, 0] state. Exploring the latter state and its stability for 0 < 5 < 0.5 (or, in fact, for any 5 since the solution 
persists V5), the corresponding Jacobian of Eq. © assumes the form 

(—1 + 8 1 —5 \ 

J = 1 —2 1 . (10) 

V -5 1 -1 + 8 ) 

The associated stability eigenvalues become Ai = {0,0}, \ 2 = + \f0ei, and A 3 = ± \/2(l — 2 8)ei. Of particular note 
is the eigenvalue pair A 3 which remains imaginary for 0 < 8 < 0.5. Thus, this solution is stable for 0 < 8 < 0.5, but 
it destabilizes through a pitchfork bifurcation, giving rise to the discrete vortex (inheriting its stability) for 8 > 0.5. 

(2) Given our analysis of the [0, 7 r, 0] state above in connection with the vortex bifurcation, we will not discuss 
further here the case with anisotropy between sites with phase 0 . 

For the [0, 7 r, 0] state with anisotropy between the site with phase n and one of the 0 phase sites, the Jacobian 
matrix is 

/0 1 - 1 \ 

J = 1 —1 — 5 5 (11) 

\-l 5 1-5/ 


from which we find the stability eigenvalues to be Ai = {0,0} and X 2 = ±y 2(—5 + V 8 2 + 3)e, A 3 = 

± \J 2(5 + V5 2 + 3 )ei. The eigenvalue pair A 2 remains real valued throughout the change of 5 and thus we see that 
this solution remains unstable. 

(3) In the isotropic lattice the stability eigenvalues of the hexagonal [0, 0,0] configuration are found to be Ai = {0, 0}, 
A 2 = ±V 6 e, and A 3 = ±\/ 6 e. Not surprisingly, due to the adjacent in-phase sites, this is unstable. In the case of 5 ^ 1, 
this instability persists. From our above analysis, the relevant eigenvalues can be directly found to be Ai = {0,0}, 
A 2 = ±\/ 6 e, and A 3 = ±\/2(25 + l)e. Despite the significant dependence of A 2 on 5, we see that this configuration is 
indeed generically expected to be unstable. 


B. Honeycomb Lattice 

We now analytically examine the stability and transformation of solutions in the anisotropic honeycomb lattice and 
its six-site cell (|M| = 6 ). Similar to what was done for the hexagonal lattice, we will deduce relationships between the 
phases for a few prototypical configurations in the honeycomb lattice. In each solution, the anisotropy lies between the 
sites with phases 63 and 64 and also between sites with phases 6 q and 63 . The structures we will present here assume 
the following form in the isotropic limit of 5=1: (1) [0, yy, 7 r, (charge 1 vortex), (2) [0, ^-,2it, Ayt] 

(charge 2 vortex), (3) [0, n, 0, n, 0, 7 r] (4) [0,0, 0,0, 0,0]. While additional configurations are possible here (in particular 
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all combinations of 0 and 7r phase are possible within the 6 sites), these configurations are the most interesting ones 
and will provide a basic understanding of the stability properties of the anisotropic system. 

(1) The honeycomb six-site charge 1 vortex (9\ = 0, |A0| = {}) satisfies the following phase relationships which can 
be deduced from Eq. 


0 3 = 

9 1+2 arccos 

92 = 

9\ + 9 3 

2 

9 4 = 

9\ + 7T 

e 5 = 

9i + 9 3 

2 + 7r 

0 6 = 

9 3 + 7T 



0.5 < 5 < 1 


( 12 ) 


(mod 27 t). 


Using these results with the Jacobian matrix of Eq. we find that the stability eigenvalues of the charge 1 vortex 
are Ai = {0,0}, A 2 = ±i/f, A 3 = A 4 = ±\/t', and 

A 5i6 = ±^j ± \J- 4+ ^ +9cos (a(5)(3(5)) +4 S 2 cos(4 a(6)(3(5)) + ( 16,5 2 ~ 9 ) cos 2 (a(J)) - 8 S 2 cos 4 (a(J))^ e 

where a(5) = 9 1 + arccos(^j), (3(6 ) = arccos (^y). It is clear that at least the eigenvalue pairs Ai, A 2 , and A 3 
all remain real-valued as <5 —> 0.5 and hence the charge 1 vortex is unstable on the entire interval 0.5 < 5 < 1. 

The changes in the relative phases (mod 27 t) as a function of <5 for 0.5 < 5 < 1 are given by: 

\0 2 - 9 4 \ = 1 0 3 - 0 2 \ = \0 5 - 9 4 \ = l^e - 0h\ = arccos 

\0 4 — 9 3 \ = tt — 2 arccos 
\9\ — 9 e \ = n + 2 arccos 



In the isotropic case (<5 = 1), all the relative phases satisfy |Ad| = J. But as 5 —>■ 0.5, we find that \9i +4 — —>• 0 

for i = 1,2,4, 5, while \9 4 — 0 3 | —>• 7r and, \9 4 — 9 6 \ —»• 7r. Also as J —)• 0.5, from Eq. (fl2l) we theoretically predict 
that the charge 1 vortex in the isotropic lattice merges at with \9 4 ,9 4 ,9 4 ,9 4 + tt,9 4 + n,9i + 7r] at 6 = 0.5. As 
was shown in more detail with the hexagonal three-site charge 1 vortex, the stability eigenvalues for the interval 
0 < 5 < 0.5 can be obtained and appear more simply as Ai = {0,0}, A 2 = ±v2e, A 3 = ±\/2(l — 25)e, A 4 = ±\Z 6 s 

and, As^ = ± (3 — 26 ± \/~45 2 + 4 5 +~9)e. Thus, we have instability throughout the change in the anisotropy. More 
generally, the established “rule of thumb” for self-focusing nonlinearities is that whenever sites of the same phase are 
adjacent to each other, the configuration will inherit an instability associated with a real eigenvalue pair. Hence, it 
is natural to expect, based on the structural form of the configuration [9 4 l 9 4 , 9 4 , 9\ + 7r, 9\ + 7r, 9 4 + 7r], that it will be 
unstable for all values of <5. Nevertheless, one of its real pairs for 5 < 0.5 will become imaginary for 5 > 0.5, giving 
rise to an unstable daughter state, namely the single charge vortex solution. 

(2) The honeycomb six-site charge-2 vortex (0 4 = 0, A0 = ^ for 5 = 1) satisfies 


9 3 = 


92 = 

9 4 = 
0 5 = 
e 6 = 


9 1+2 arccos ( — 

9\ + 9 3 
2 

6»i +2 t r 

9\ + d 3 


2 

%+2tt 


■ 2n 



0.5 < 5 < 1 


(13) 


(mod 47 t) 


and the stability eigenvalues for the charge 2 vortex are Ai = {0,0}, \ 2 = A 3 = A 4 = ±and 

A5.6 = A t 


{ 


\ 


1 + 2 S S2 ± J —4 + + 9 cos(a(— 5)(3(—5)) + 46 2 cos(4a(— 5)/3(—5)) 


' 16(5 2 -9 


) cos 2 (a(— 5)) — 8<5 2 cos 4 (a(— 5)) I ei, 
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where as before a(<5) = 0\ + arccos (^j), /3(d) = arccos (^j)- 

The changes in the relative phases as a function of <5 for 0.5 < <5 < 1 are given by: 

|02 ~ 0i| = |0 3 - d 2 \ = \d 5 - 0 4 1 = |06 - 0 5 1 = arccos 

104 — 0 3 | = 2 n — 2 arccos 

|0i — 06| = 27T + 2 arccos . 

In the isotropic honeycomb lattice (6 = 1), all the relative phases satisfy |A0| = But as <5 —0.5, we find that 
|0;+i — 6 i\ —>■ 7T for i = 1,2,4, 5, while |04 — 0 3 | —)• 0 and, |0i — 0§\ -A 0. Also, as S —>■ 0.5, from Eq. (fl3l) we 
theoretically predict that the charge 2 vortex in the isotropic lattice merges at <5 = 0.5 with [0i,0i + 7r,0i + 27 t,0i + 
27t,0i + 37T,0i + 47 t] = [0i,0i + 7r,0i,0i,0i + 7T,0i] (mod 27r). The stability eigenvalues of the latter state for 

0 < 6 < 0.5 are Ai = {0,0}, A 2 = ±y/2ei, A 3 = ±^2(1 - 2S)ei, A 4 = ±V 6 si, A 5 = ±^/(-3 + 25 - V4 S 2 + 46~+9)£ = 

± \J (3 — 2d + y/Z5 2 + 46 +"9 )ei, A6 = ± \J (—3 + 26 + \/4<5 2 + 46 + 9)e. Note that A6 is real for 0 < 6 < 0.5 and so 
it is expected that the stable vortex of the isotropic limit with <5=1 becomes unstable for some <5 in 0.5 < <5 < 1. 
As regards the state [0i,0i + 7r,0i,0i,0i + 7r,0i] (mod 27r), the above analysis predicts that it also undergoes a 
supercritical pitchfork bifurcation giving rise to the charge 2 vortex which inherits its (in)stability properties and is 
eventually full stabilized for a larger value of <5 in 0.5 < 6 < 1. In fact, using the expression for A6 given above for 
the charge-2 vortex, we find that the relevant critical point is 6 = 0.716, which we will compare with our numerical 
computations in the following section. 

(3) In a similar way to the analysis done above, the stability eigenvalues for the six-site configuration of al¬ 
ternating sites [ 0 , 7 r, 0 , 7 r, 0 , 7 r] are found to be Ai = { 0 , 0 }, \ 2 = ±\/ 6 ei, A 3 = ±\/ 2 ei, A 4 = ±\/2(2 6 + 1 )ei 

A 5_6 = ± \J(3 + 26 ± \/4A 2 — 46 +"9 )ei. The imaginary eigenvalues indicate that this configuration is stable both 
in the isotropic and anisotropic lattices. 

(4) Finally the six-site configuration of in-phase sites [0,0, 0,0, 0,0] has the stability eigenvalues Ai = {0,0}, X 2 = 

±\/ 6 e, A 3 = ±\/2e, A 4 = ± \/2(26 + l)e As^ = ± \J (3 + 25 ± \/4 5 2 — 46 + 9)e. Given the presence of real eigenvalues, 
this configuration is unstable. Similarly to what we saw previously for the 3-site configuration in the hexagonal lattice, 
the presence of adjacent in-phase sites is detrimental to the stability of this configuration for arbitrary values of 6 . 
Hence, no stabilization of the relevant state is anticipated, independently of the particular value of the anisotropy 
parameter 6 . 


III. NUMERICAL RESULTS 


In this section we present our numerical findings for the various configurations in the hexagonal and honeycomb 
lattices and compare these results with the theoretical findings from the previous section. In all cases, we use a 
Newton-Raphson fixed point iteration to identify the full numerical solution over the two-dimensional lattices. This 
process is initiated at the AC limit and continued to a small coupling while maintaining <5 = 1 to yield the isotropic 
case. In most of the examples we consider, the Newton-Raphson interation is continued to e = 0.01, but a few cases 
will also be examined where the fixed point iteration is continued to a higher value of e. The anisotropy (as described 
in a previous section) is introduced by letting <5 deviate from the isotropic unity value and performing a continuation 
in decreasing values of the parameter towards <5—^0. We present figures that show each configuration, along with its 
phase portrait and spectral plane at some 6 before and after the relevant bifurcation points, comparing the latter with 
our theoretical predictions. It is generally found that for values of the coupling, e, near the AC limit the theoretical 
predictions match the numerical results very well. 

The theoretical predictions of the linearization eigenvalues will be compared to the numerical results for the linear 
stability of the stationary solution v m ^ n exp(iz) by using the ansatz 


■7(0 


_A z 


b* 


.A 2) 


(14) 


The eigenvalue problem that follows is then solved for the eigenvalues A and eigenvectors (a m , n , b m ,n) T ■ The asterisk 
denotes complex conjugate while T denotes transpose. 
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A. Hexagonal Lattice 

(1) We start with the hexagonal lattice and begin by presenting the results of the continuation for the three site 
charge 1 vortex with coupling e = 0.01 (see Fig. [[]). Initially, at 8 = 1, the configuration has phases 0 3 = 0,02 = 
27 t/3, 03 = 47t/ 3, with the anisotropy to be activated between the sites with phases 6 \ and 0 3 . The top row of Fig. |T| 
shows the square modulus of the field (|un,m| 2 ) for the vortex at <5 = 0.8 (left panel) and at S = 0 (right panel), while 
the second and third rows show, respectively, the corresponding phase portraits and spectral planes. The left panel of 
the fourth row shows the change in the relative phases (|A0|) and the right panel traces the linear stability eigenvalues 
with the change in anisotropy. In both images the theory (dash-dot lines) compares extremely well with the numerical 
(solid lines) results, over the entire interval of continuation of the anisotropy parameter 8. In the case of <5 —^ 0 it 
is evident that the configuration has changed its character into a [0, n, 0] configuration, as theoretically predicted. 
The relative phase and eigenvalue predictions of the bottom row indeed confirm that the pitchfork bifurcation takes 
place at 8 = 0.5 and leads to a collision with the configuration [—7t/3, 27 t/ 3, — 7r/3] at 8 = 0.5 (which as theoretically 
predicted has the form [0i, 0 3 +7r, 0 i], i.e., up to a trivial phase is a [0,7r, 0] configuration). We note in passing that we 
have also examined the relevant continuation and bifurcation for other values of £ (such as e = 0.05), finding similar 
qualitative results (although quantitative details, such as the critical value of 8 do change). 

(2) The hexagonal [0,7r, 0] configuration can be seen in Figs. [2] and [3] In the isotropic lattice of e = 0.01 this 
solution is unstable. However stability can be achieved if the anisotropy is activated between the the two nodes with 
phase 0, while any other placement maintains the instability. This is because if the nodes with the same phase are 
connected at the 8 —> 0 limit (where we can think of the three nodes as effectively being on a straight line), then, 
as discussed above, the instability due to a real eigenvalue pair will be maintained. In the case where a stabilization 
effect is observed (in Fig. [2]) [0,7r, 0], as 8 decreases, a weakened bond between the 0 phase nodes results and eventually 
an effective [0,7r, 0] state along a line is effectively obtained which is well known to be stable, as a one-dimensional 
configuration, for small s (TV. The change of stability occurs at 8 = 0.49, in excellent agreement with the theoretical 
prediction of 8 = 0.5. On the other hand, an anisotropic weakening of the bond between the site with phase 7r and 
either of the 0 phase sites does not bring about stability and eventually just yields an unstable waveform (along an 
effective line) [0,0, tt], which has been demonstrated to be unstable in Id settings (l9| . 

(3) The hexagonal [0,0,0] configuration (Fig.Qj) is, not surprisingly, unstable in the isotropic lattice (e = 0.01) due 
to the adjacent in-phase sites and, in full accordance with our theoretical predictions, this remains unstable as <5 —> 0. 


B. Honeycomb Lattice 

We now discuss the effects of anisotropy in the honeycomb lattice. For the case of the six-site charge 1 vortex (see 
Fig. [5j), we again see a very good comparison between the numerical results and the theory. The unstable charge 1 
vortex remains unstable throughout the anisotropic variation of 5 in the interval [0,1). At about 8 = 0.70 one of the 
real eigenvalue pairs As^ collides with the origin of the spectral plane and becomes purely imaginary. At 8 = 0.5 the 

pair which was theoretically predicted to be A 3 = ± \J ( 4| ^ A]£ collides with the pair Ai, giving rise to the bifurcation 
that was theoretically predicted to arise at this critical point. Indeed this bifurcation transforms the vortex for 8 < 0.5 
into the unstable [7 t/3, 7t/3, 7t/3, 27t/3, 2n/3, 27t/ 3] honeycomb configuration (i.e., a [0, 0,0,7r, 7r, 7r] state up to a trivial 
phase shift, as discussed in section II). 

The stable honeycomb charge 2 vortex (shown in Fig. [6]) becomes unstable at approximately 5 = 0.70, in very good 
agreement with the theoretical prediction of 8 = 0.716; cf. the relevant discussion in section II. The instability arises 
due to a pair of eigenvalues from As^ becoming real-valued. Subsequently, as theoretically predicted, for 8 = 0.5, 
a pitchfork bifurcation eventually transforms the vortex into the state [— 7 t/3, 27t/3, — 7t/3, —7t/3, 27t/ 3, —7r/3], (i.e., a 
[0,7T, 0,0,7r, 0] state up to a trivial phase shift). The latter configuration has the same stability characteristics for 
8 < 0.5 that the vortex state possesses in the vicinity of 8 > 0.5. Therefore, it possesses a single real eigenvalue pair 
as confirmed in the right panels of Fig. [6] Naturally, for 8 > 0.5, this [0,7r, 0, 0,7r, 0] persists but acquires a second 
real pair as a result of the supercritical pitchfork bifurcation. 

In addition to these vortex configurations, in Figs. [7] and [5J we also examined the states with phase configurations: 
[0,7r, 0,7T, 0,7r] and [0,0, 0,0, 0,0], respectively. The former, as expected (for this small e and given its alternating 
phase structure) is found to be linearly stable for all the considered values of 8 , while the latter is found to be highly 
unstable bearing 5 distinct real eigenvalue pairs (as anticipated due to the presence of sites of the same phase adjacent 
to each other). 




IV. DYNAMICS 


In this section we numerically examine the nonlinear dynamics of the unstable solutions discussed in the previous 
sections. Each unstable solution is perturbed slightly in the direction of the eigenvector corresponding to the most 
unstable eigenvalue, in order to seed the relevant instabilities. A fourth order (explicit) Runge-Kutta algorithm (RK4) 
has been used in order to obtain the relevant dynamical evolution results. It is observed that the coupling significantly 
controls the nature of the dynamical evolution. This is natural since a larger coupling allows nearest neighbors to 
interact more strongly. In the hexagonal lattice, we show the evolution at a fixed background coupling of e = 0.01 
and also e = 0.2. In all cases, the smaller coupling leads to a robust (multi-site) breather form, while the larger 
coupling produces a single robust site. For the honeycomb lattice case, a coupling of e = 0.2 is used. All the unstable 
honeycomb six-site solutions evolve into multi-site breathers (where the number of sites participating with a large 
norm in the final configuration varies from case to case; see below). In either lattice we hold the anisotropy fixed at 
8 = 0.8 for instabilities above the critical threshold of 8 = 0.5, and to 8 = 0.2 if the instability occurs below 8 = 0.5. 
The one exception is for the charge 2 vortex where we present the dynamics at <5 = 0.6 (instead of 8 = 0.8), i.e. after 
the onset of instability but before the critical transformation point of that state. 


A. Hexagonal Lattice 

The solutions found to be unstable in the hexagonal geometry are (1) [0,7r, 0] with the anisotropy between the 0 
phase sites when 8 > 0.5 (due to the bifurcation of the vortex state), (2) [0, tt, 0] with the anisotropy between the site 
with phase n and one of the sites with phase 0, (3) [0, 0,0]. 

(1) Figures 151 and HU1 exhibit the dynamics in the anisotropic hexagonal lattice for [0,7r, 0] when the anisotropy is 
prescribed to be between the sites with phase 0 with e = 0.01 and e = 0.2 respectively. In both instances, we use 
6 = 0.8. In the case of the weaker coupling, the propagating solution shows oscillatory (i.e., breathing) behavior. In 
fact, we observe similar features for the weak coupling case of e = 0.01 for all the unstable solutions studied here. For 
the larger coupling of e = 0.2 in Fig. [10] we see a clear destruction of the original waveform and an emergence of a 
single surviving site that persists. It may be interesting in future work to explore the transition regime between weak 
and strong coupling and the associated implications for the nature of the resulting states. 

(2) The [0,7r, 0] solution with anisotropy between the site with phase 7r and one with 0 is unstable for all values 
of 8, as discussed previously. As in the previous example, for e = 0.01 the dynamics shows an oscillatory movement 
during propagation. However, when the coupling is set to e = 0.2, destruction of the wave is observed with, as before, 
a single site persisting for long times. Given the similarity of these findings to those of Figs. [9l and [TUI we omit them 
here. 

(3) Finally the evolution of the form [0,0, 0] is seen in Figs. |TT|and[T2] The solution is unstable for all 0 < 8 < 1 . 
We see essentially the same qualitative behavior as for the previous two cases, as regards the asymptotic fate of the 
unstable waveforms. We have also checked that this phenomenology arises for different values of 8 , such as 8 = 0.2 
(results not shown here). 


B. Honeycomb Lattice 

Finally, we now turn to case examples of the dynamical evolution on the honeycomb lattice. The unstable solutions 
that we study in this case are: (1) the charge one vortex, (2) the charge-2 vortex for 8 < 0.716 (below 0.5, recall that 
this states morphs into an unstable [0, 0 + n, 0,0, 0 + n, 0] state), (3) the in-phase solution [0, 0,0, 0,0, 0]. Given that 
the results for weak coupling are similarly (breathing) as in the previous subsection, we focus on the case of the larger 
coupling e = 0.2. 

(1) The charge 1 vortex is unstable for the full range of the anisotropy considered herein. In Fig. [T3] we explore 
its unstable dynamics for a coupling strength of e = 0.2 and 8 = 0.8; similar results have also been found for other 
values of <5, as e.g. in Fig. [14] for 8 = 0.2. The dynamics yields a multi-site excitation, but with a repartitioning of 
the relevant intensity so that some sites are dominant in amplitude in comparison to others. 

(2) In Fig. [15] we see the dynamics of the charge 2 vortex after the onset of instability, at <5 = 0.6. The result is a 
six site breathing structure with a complex norm redistribution. In Fig. 1161 we show the dynamics of the unstable 
solution of the form [0,0 + 7r, 0, 0,0 + 7r, 0], resulting from the pitchfork bifurcation of the charge-2 vortex for the case 
of 8 < 0.5. Specifically, in this example for 8 = 0.2, a breathing six-site excitation appears to persist. 

(3) Finally, the dynamical evolution of the unstable state [0, 0,0, 0,0, 0] for 8 = 0.8 in Fig. [T7] (but also for other 
values of 8) illustrates the dynamical tendency of this state towards configurations with fewer -arguably two, at the 
final evolution snapshot shown- dominant (in amplitude) sites. 
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V. CONCLUSIONS AND FUTURE CHALLENGES 

In summary, in the present work, we have explored the existence, stability and dynamics of localized states (focusing 
on multi-site solitonic and vortex states) in hexagonal and honeycomb lattices. We considered the prototypical unit 
cells in each case, namely a 3-site one in the hexagonal case and a 6-site one in the honeycomb case. Analytical con¬ 
siderations in the vicinity of the anti-continuum limit permitted us to identify the states in the presence of anisotropy 
in an approximate analytical form and gave us the ability to consider the linear (spectral) stability eigenvalues and 
obtain approximate analytical expressions for them. These results allowed us to elucidate the pitchfork bifurcations 
that lead to the disappearance of states such as vortices (and the destabilization of other solitonic states) as the path 
from more effectively one-dimensional to effectively two-dimensional configurations is traversed. These existence and 
stability findings were also found to be in good agreement with detailed numerical continuations (over the anisotropy 
parameter), at least for small values of the coupling. Finally, the dynamics of the relevant structures were examined, 
allowing us to identify some gross features, including the breathing nature of the instability for very weak e and the 
potential for stronger localization (typically to a smaller number of sites) ensuing as a result of instability for stronger 
e. 

A significant number of possibilities emerge from the present work for future explorations. On the one hand, it would 
be relevant and interesting to examine in more detail the dynamical evolution scenarios of the model, and to provide 
a more systematic characterization of the propagation outcomes for cases of both weaker and stronger coupling. On 
the other hand, extending similar studies to the case of Kagome lattices and their flat bands, identifying the spectral 
properties not only of the solitons/vortices [lj| but also of the compactly supported structures identified therein 
would be a timely theme. Finally, extending such considerations to three-dimensional lattices of different types would 
also pose significant new challenges and can be expected to feature intriguing bifurcation phenomena and states of 
interest. Efforts along these directions are presently underway and will be reported in future publications. 


10 





1 

1 

i 

1 

0.5 


0.5 



8 


O 

0 

<-T 0 

o 

oT 0 

0 


8 


0 

-0.5 


-0.5 


-1 

i 

-1 

1 


-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 

X X 




FIG. 1: Hexagonal three-site charge 1 vortex i.e. the hexagonal [0, 27t/3, 47t/3] configuration (at e = 0.01). The anisotropy 
here is activated between the sites with (initial) phases 0 and 47r/3 when 5=1. The top row displays the modulus squared of 
the configuration corresponding to anisotropic parameter 5 = 0.8 (left panel) and 5 = 0 (right panel). The second row shows 
the phase portraits and the third row shows the spectral plane for the same values of the anisotropy, 5 = 0.8 (left), and 5 = 0 
(right). In the fourth row, the left panel shows the comparison between the theoretical (dash-dot lines) and numerical (solid 
lines) changes in the relative phases. The charge 1 vortex collides at 5 = 0.5 with the stable (for lower values of 5) hexagonal 
[— 7t/3, 27t/3, — 7t/3] configuration. When 5 = 0 this is equivalent to the [—7r/3, 27t/3, —7r/3] configuration, i.e., effectively a 
[0, 7r, 0] configuration. The fourth row right panel shows a theoretical (dash-dot lines) versus numerical (solid lines) comparison 
of the stability eigenvalues for 0 < 5 < 1. 
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FIG. 2: Hexagonal three-site [0,7r, 0] configuration (e = 0.01). The top row displays the modulus squared of the field at 5 = 0.8 
(left column) and 5 = 0 (right column). The second row shows the corresponding phase portraits while the third row displays 
the respective spectral planes. In this case, the anisotropy is invoked between the two nodes with phase 0. Hence, as discussed 
in the text a stabilization is observed for 5 < 0.5. 
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FIG. 3: Hexagonal three-site [0, n, 0] configuration (e = 0.01). The top row displays the modulus squared of the field with 
8 = 0.8 (left column) and 6 = 0 (right column). The anisotropy is invoked between the node with phase 7r and one of the 
nodes with phase 0. At <5 = 0 the result is equivalent to an effective [0,0, tv] configuration, along a line. Hence, in this case the 
instability is preserved for all values of 8. 
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FIG. 4: Hexagonal three-site [0,0,0] configuration (e = 0.01). The top row displays the modulus squared of the field with 
anisotropy 0.8 (left column) and 0 (right column). The configuration is found to be unstable for all values of 8, as also predicted 
theoretically. 
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FIG. 5: Honeycomb six-site charge 1 vortex i.e. the honeycomb [0, 7t/3, 27t/3, n, 47t/3, 57t/3] configuration at the isotropic limit 
(e = 0.01). The top row displays the modulus squared of the configuration corresponding to anisotropic parameter 5 = 0.8 
(left panel) and 5 = 0.3 (right panel). The second row shows the phase portraits and the third row shows the spectral plane 
for the same values of the anisotropy, 5 = 0.8 (left), and S = 0.3 (right). In the last row the left panel shows the comparison 
between the theoretical (dash-dot lines) and numerical (solid lines) changes in the relative phases. The vortex collides at 5 = 0.5 
with the [7r/3, 7r/3, 7r/3, —27r/3, —27r/3, —27r/3] configuration. At 8 = 0, for the present form of anisotropy, this is equivalent 
to the two configurations along a line, namely [7t/3, vr/3, vr/3] and [— 27t/3, — 27t/3, — 27r/3]. The bottom right panel shows the 
comparison of the theoretical (dash-dot lines) versus numerical (solid lines) linear stability eigenvalues for 0 < S < 1. 

















14 



1 

0.5 
-T 0 
-0.5 
-1 

-0.01 -0.005 0 0.005 0.01 

X 





1 

0.5 

0 


-0.5 


-1 




0 

8 

o o o 

8 
0 




-0.1 -0.05 


0 


0.05 0.1 




FIG. 6: Honeycomb six-site charge 2 vortex i.e. the honeycomb [0, 27r/3, 47 t/3, 27t, 87t/3, 10-7t/ 3] configuration (e = 0.01). The 
top row displays the modulus squared of the configuration corresponding to anisotropic parameter 8 = 0.8 (left panel) and 
8 = 0.3 (right panel). The second row shows the phase portraits and the third row shows the spectral plane for the same values 
of the anisotropy, 8 = 0.8 (left), and <5 = 0.3 (right). In the last row the left panel shows the change in the relative phases. The 
charge 2 vortex collides at 8 = 0.5 with the [—tt/3, 27t/3, —7t/3, —7t/3, 2tt/3, — 7t/ 3] configuration. The bottom right panel shows 
the theoretical (dash-dot lines) versus numerical (solid lines) comparisons of the linear stability eigenvalues for 0 < 8 < 1. 
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FIG. 7: Honeycomb six-site out-of-phase configuration, i.e. honeycomb [0, n, 0, n, 0,7r] configuration (e = 0.01). The top row 
displays the modulus squared of the field with <5 = 0.8 (left panel) and 8 = 0.3 (right panel). At <5 = 0 this is equivalent to the 
two configurations along a line of the form: [0, tv, 0] and [n, 0,7r]. Hence it retains its stability throughout the interval 0 < 5 < 1. 
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FIG. 8: Honeycomb six-site in-phase configuration, i.e. honeycomb [0, 0, 0, 0, 0, 0] configuration (e = 0.01) The top row displays 
the modulus squared of the field with anisotropy 0.8 (left column) and 0.3 (right column). As expected by the adjacency of 
sites with the same phase, the anisotropy cannot prevent this configuration from being highly unstable. 
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